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Quantum  chemistry-based  dipole  polarizable  and  nonpolarizable  force  fields  have  been  developed 
for  1, 3, 5-triamino-2, 4, 6-trinitrobenzene  •TATB*.  Molecular  dynamics  simulations  of  TATB  crystals 
were  performed  for  hydrostatic  pressures  up  to  10  GPa  at  300  K  and  for  temperatures  between  200 
and  400  K  at  atmospheric  pressure.  The  predicted  heat  of  sublimation  and  room-temperature 
volumetric  hydrostatic  compression  curve  were  found  to  be  in  good  agreement  with  available 
experimental  data.  The  hydrostatic  compression  curves  for  individual  unit  cell  parameters  were 
found  to  be  in  reasonable  agreement  with  those  data.  The  pressure-  and  temperature-dependent 
second-order  isothermal  elastic  tensor  was  determined  for  temperatures  between  200  and  400  K  at 
normal  pressure  and  for  pressures  up  to  10  GPa  on  the  300  K  isotherm.  Simulations  indicate 
considerable  anisotropy  in  the  mechanical  response,  with  modest  softening  and  significant  stiffening 
of  the  crystal  with  increased  temperature  and  pressure,  respectively.  For  most  properties  the 
polarizable  potential  Was  found  to  yield  better  agreement  with  available  experimental  properties. 
©  2009  American  Institute  of  Physics.  *doi:10.  1063/1 .3264972* 

I.  INTRODUCTION 

1, 3, 5-triamino-2, 4, 6-trinitrobenzene  ♦TATB*  is  an  ener¬ 
getic  material  used  in  high-performance  military  applica¬ 
tions.  The  molecular  structure  of  TATB,  shown  in  Fig.  1,  is 
characterized  by  strong  intramolecular  hydrogen  bonding  be¬ 
tween  neighboring  amine  and  nitro  groups.  TATB  exhibits 
remarkably  low  sensitivity  to  accidental  initiation  by  unin¬ 
tended  external  stimuli  such  as  shock,  thermal,  or  electrical 
loading.1,2  Numerous  spectroscopic,  chemical,  and  structural 
studies  have  been  performed  in  attempts  to  understand  the 
origin  of  this  stability.  Among  the  explanations  that  have 
been  suggested  and  explored  are  chemical  stabilization  due 
to  strong  intramolecular  push-pull  conjugation3  and  related 
electronic  interactions;4  stabilization  due  to  intermolecular 
hydrogen  bonding;2  centrosymmetry,  or  lack  thereof  based 
on  second-harmonic  frequency -doubling  experiments;5*8  and 
the  low  potential  energy  barrier  to  compression  and  glide 
between  layers  in  the  TATB  crystal  structure 9 

Cady  and  Larson10  determined  that  at  room  temperature 
and  atmospheric  pressure  TATB  crystallizes  in  a  triclinic  unit 
cell  containing  two  molecules  in  the  PI  space  group,  The 
crystal  consists  of  graphitic-like  sheets  in  the  a-b  crystal 
plane,  with  significant  inter  and  intramolecular  hydrogen 
bonding  in  that  plane  *see  Fig.  1*.  By  contrast,  interplanar 
interactions  along  the  c-axis  are  limited  to  weak  van  der 
Waals  forces.  These  distinct  differences  in  bonding  along 


different  crystal  directions  lead  to  significant  thermome¬ 
chanical  anisotropy,  and  numerous  experimental  and  theoret¬ 
ical  studies  have  been  performed  to  characterize  and  explain 
the  effects  of  this  anisotropy  on  observed  TATB  properties. 
The  anisotropic  thermal  expansion  of  TATB  was  studied  by 
Rizzo  and  Kolb;1  those  workers  reported  an  order-of- 
magnitude  difference  between  the  linear  expansion  coeffi¬ 
cients  along  a  and  b  compared  to  c.  Olinger  and  Cady12 
conducted  an  experimental  investigation  of  room- 
temperature  TATB  crystal  structure  under  hydrostatic  com¬ 
pression  from  normal  pressure  to  7.47  GPa.  In  order  to  ob- 


•"Electronic  mail:  d.bedrov@uuh.edu.  FIG.  1.  Molecular  and  crystal  structure  of  TATB. 
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tain  the  lattice  parameters  from  their  data,  0 linger  and  Cady 
assumed  that  the  ratio  a/b  was  constant  and  also  that  points 
in  adjacent  a-b  planes  maintained  fixed  relative  positions 
under  compression.  Recently,  Stevens  et  al.13  reported  volu¬ 
metric  data  for  hydrostatic  compression  of  TATB  crystal  for 
pressures  up  to  13  GPa.  Their  results  were  generally  consis¬ 
tent  with  those  of  Olinger  et  al.,12  although  Stevens  et  al.13 
did  report  a  slightly  stiffer  mechanical  response.  Interest¬ 
ingly,  Stevens  et  al.  1  observed  a  small  kink  in  the  isotherm 
at  P  =  8  GPa  when  plotting  the  data  in  the  pseudo  shock- 
velocity/particle-velocity  plane.  Pravica  and  co-workers14,15 
studied  the  infrared  spectrum  of  TATB  under  compression  to 
40  GPa  in  a  diamond  anvil  cell  using  Fourier  Transform 
Infrared  Spectroscopy.  They  investigated  shifts  in  spectral 
peaks  in  the  far-  and  mid-infrared  regions  and  observed  de¬ 
creases  in  N-H  and  N-0  stretching  frequencies  with  increas¬ 
ing  pressure,  indicative  of  increased  hydrogen  bonding  under 
compression.  They  obtained  data  during  both  loading  and 
unloading  and  found  no  evidence  for  either  phase  transitions 
or  chemical  decomposition  in  the  pressure  interval  studied. 
Fits  of  experimental  volumetric  compression  data  to  various 
equations  of  state  *EOSs*  have  provided  estimates  of  isother¬ 
mal  bulk  modulus  *K0*  and  its  initial  pressure  derivative 
•K~.  However,  to  the  best  of  our  knowledge  no  attempts 
to  resolve  the  complete  set  of  elastic  constants  for  TATB 
crystal  have  been  made  from  either  experiment  or  theory. 

Theoretical  and  simulation  studies  of  TATB  compression 
have  also  been  reported.  Pastine  and  Bernecker16  reported  a 
theoretical  EOS  obtained  from  calculations  based  on  simple 
Lennard-Jones  interactions  that  only  addressed  interplane  in¬ 
teractions  along  the  c-axis.  This  EOS  predicted  a  compress¬ 
ibility  significantly  lower  than  values  based  on  experiment  3 
Byrd  and  Rice17  used  plane  wave  density  functional  theory 
•OFT*  calculations  to  investigate  TATB  crystal  under  hydro¬ 
static  compression.  This  work  showed  that  OFT  predictions 
for  energetic  crystals  are  quite  sensitive  to  the  size  of  basis 
set  and  that  DFT  provides  a  poor  description  of  unit  cell 
lattice  vectors  at  lower  pressures,  presumably  due  to  neglect 
of  dispersion  interactions  in  the  theoretical  model. 

Sewell  reported  studies  of  temperature-dependent  unit 
cell  volume,  lattice  parameters,  and  energies  obtained  from 
isobaric-isothermal  *NPT*  Monte  Carlo  *MC*  simulations  of 
TATB  crystal,  based  on  a  simple  Buckingham-plus-charges 
intermolecular  interaction  potential  and  extension  to  treat 
atom-centered  multipoles  through  octupole.  The  effects  of 
system  size,  intramolecular  flexibility  -fully  rigid  molecules 
versus  ones  that  included  hindered  N02  rotations*,  and  order 
of  expansion  of  atom-centered  multipoles  were  considered. 
Although  inclusion  of  atom-centered  multipoles  yielded  im¬ 
proved  fits  to  electrostatic  potentials  compared  to  results  ob¬ 
tained  using  simple  partial  atomic  charges,  they  did  not  lead 
to  significantly  different  temperature  dependencies  for  en¬ 
ergy  or  density.  By  contrast,  inclusion  of  N02  rotation  was 
found  to  have  a  significant  effect  on  both  thermal  expansion 
and  specific  heat. 

Gee  et  al.  1  developed  a  fully  atomistic  flexible  molecule 
force  field  *FF*  and  used  it  in  molecular  dynamics  *MD* 
simulations  to  study  pressure-volume-temperature  properties 
of  TATB  crystal.  In  their  FF  Gee  et  al.1'  obtained  Lennard- 
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Jones  nonbonded  parameters  by  fitting  to  gas-phase  TATB 
dimer  energies  obtained  from  quantum  chemistry  *QC*  cal¬ 
culations  using  the  MP2/6-31G*d,p*  theory  and  basis  set 
combination  for  geometries  relevant  to  the  crystal.  Intramo¬ 
lecular  covalent  bond  and  three-center  angle-bending  inter¬ 
actions  were  taken  from  the  DREIDING  FF  while  torsional 
potentials  for  N02  and  NH2  groups  were  parameterized  to 
reproduce  previously  published  QC  calculations  of  the  dihe¬ 
dral  barrier  heights.  Gee  et  al.19  reported  identical  values  of 
the  Lennard-Jones  well  depth  parameters  for  C-C  and  N-N 
nonbonded  interactions.  The  optimized  value,  0.01  kcal/mol, 
is  six  times  lower  than  the  well  depth  they  used  for  H-H 
interactions  and  is  in  striking  disagreement  with  correspond¬ 
ing  values  used  in  other  organic  molecule  FFs  such  as 
OPLS-AA20  or  AMBER.21  Nevertheless,  Gee  et  al.19  pre¬ 
dicted  crystal  structure,  temperature-  and  pressure-dependent 
lattice  parameters,  and  other  structural  metrics  that  are  in 
good  agreement  with  the  experimental  data  of  Olinger  and 
Cady. 2  However,  it  was  reported  subsequently22  that  the 
heat  of  sublimation  for  the  Gee  et  al. 9  FF  is  •  H^ 
=  31.1  kcal/mol,  which  is  significantly  lower  than  reported 
experimental  values  due  to  Garza23  *43.1  kcal/mol*  and 
Rosen  and  Dickinson24  *40.2  kcal/mol*.  Taken  together, 
these  factors  suggest  that  the  Gee  et  al.19  FF  may  signifi¬ 
cantly  underestimate  the  interaction  strength  between  TATB 
molecules  in  the  crystal.  However,  another  key  result  of  the 
work  of  Gee  et  al.  is  that  MD  predictions  for  TATB  are 
highly  sensitive  to  details  of  the  FF  employed.  In  particular, 
they  showed  that,  starting  from  the  experimental  structure, 
attempts  to  simulate  TATB  crystal  using  generic  FFs  such  as 
CVFF,25  CFF91-950*,26  COMPASS,2  DREIDING,28  and 
UFF  *Ref.  29*  resulted  in  phase  transitions  into  polymorphs 
that  are  not  observed  experimentally.  Gee  et  al.  also  reported 
a  study  using  their  FF  in  which  they  investigated  adhesion 
energies  between  TATB  crystal  and  various 
fluoropolymers.30  In  this  work  standard  combining  rules 
were  used  to  obtain  nonbonded  interactions  between  TATB 
and  the  polymers.  Taking  into  account  the  likely  underesti¬ 
mation  in  the  Gee  et  al. 9  TATB  FF  of  nonbonded  parameters 
discussed  above,  we  believe  it  is  difficult  to  have  high  con¬ 
fidence  that  TATB-polymer  interactions  would  be  repre¬ 
sented  accurately  in  those  simulations. 

Recently,  Rai  et  al.22  extended  the  transferable  potentials 
for  phase  equilibrium  *TraPPE*  FF  to  model  crystalline 
TATB.  This  FF,  in  our  opinion,  uses  physically  more  reason¬ 
able  parameters  for  dispersion  interactions  and  has  a  signifi¬ 
cantly  different  distribution  of  atomic  partial  charges  than  the 
FF  of  Gee  et  al. 9  The  pressure- volume-temperature  behavior 
for  TATB  crystal  predicted  from  MD  simulations  using  this 
FF  was  in  as  good  agreement  with  experimental  data  as  that 
predicted  using  the  FF  model  of  Gee  et  al.  9  However,  the 
heat  of  sublimation  predicted  from  this  FF  is  only  25.8  kcal/ 
mol;  this  value  is  almost  40%  lower  than  the  experimental 
values,  indicating  that,  despite  predicting  accurately  the  tem¬ 
perature  and  pressure  dependencies  of  the  lattice  parameters, 
the  strength  of  intermolecular  interactions  appears  to  be  un¬ 
derestimated  in  this  FF. 

In  the  present  work  we  report  results  of  molecular  simu¬ 
lations  of  crystalline  TATB  using  FF  models  based  on  the 
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Atomistic  Polarizable  Potential  for  Liquids,  Electrolytes,  and 
Polymers  *APPLE&P*3  model.  In  contrast  to  the  Gee  et  al.'9 
and  TraPPE  FFs  which  use  fixed  partial  atomic  charges,  in 
this  work  we  used  a  polarizable  FF  that  accounts  explicitly 
for  induced  polarization  effects  that  typically  exist  in  con¬ 
densed  phase  organic  crystals.  We  show  that  simulations  us¬ 
ing  this  polarizable  FF  accurately  predict  the  structure  and 
the  heat  of  sublimation  of  TATB  crystal,  A  nonpolarizable 
version  of  the  FF  has  also  been  developed;  we  demonstrate 
that  it  yields  an  adequate  description  of  the  structural  and 
thermodynamic  properties  of  the  crystal  at  considerably  re¬ 
duced  computational  expense.  Simulation  results  for 
pressure-  and  temperature-dependent  properties  are  com¬ 
pared  to  experimental  data,  including  newly  available  results 
for  unit  cell  parameters  that  complement  previously  reported 
experimental  data  due  to  Stevens  et  al.  3  In  addition  to  those 
FF  validation  results,  we  report  complete  second-order  iso¬ 
thermal  elastic  tensors  for  a  range  of  pressures  and  tempera¬ 
tures.  To  our  knowledge,  this  is  the  first  published  report  of 
the  elastic  tensor  for  TATB. 

II.  SIMULATION  AND  FORCE  FIELD  DETAILS 
A.  Form  of  the  potential  energy  function 

In  this  work  the  potential  energy  U,ot  for  a  collection  of 
TATB  molecules  is  expressed  as: 


NB»f  •  +  •  ub0nd*ril 

•  +  •  ubend**ijk* 

4  bonds 

bends 

•  y  dihedral..  ^ 

♦  •  UlmP**  $*, 

•1* 

dihedrals 

improper 

dihedrals 

where  the  sums  for  intramolecular  terms  are  over  all  covalent 
bonds,  bends,  dihedrals,  and  improper  dihedrals  ♦out-of- 
plane  bends*  in  the  system.  These  individual  contributions  to 
the  potential  energy  are  described  by  harmonic  or  trigono¬ 
metric  functions.  These  forms  are  frequently  used  in  MD 
simulations;  their  exact  forms  can  be  found  in  our  previous 
publications.33  The  FF  form  is  redundant  in  that  there  are 
more  than  3N-6  covalent  terms  per  molecule;  in  particular, 
for  each  TATB  molecule  there  are  terms  for  24  covalent 
bonds,  36  bend  angles,  48  dihedral  angles,  and  12  out-of¬ 
plane  bends. 

The  nonbonded  energy  UNB*r*  consists  of  the  sum  of 
two-body  repulsion  and  dispersion  terms  UR0*r\  the  energy 
due  to  interactions  of  fixed  partial  atomic  charges  Ucoul*r*, 
and  the  polarization  energy  u^'t*  arising  from  the  interac¬ 
tion  of  induced  dipoles  with  fixed  charges  and  other  induced 
dipoles, 

_  yRO,^,  +  ycoul#j.,  ^  (JpoLj** 


=  •  A.,  exp**  B..r,,**  C..rjj6 


'•j 


•JI&_».o.5.  .,c 


B..rtl 


i*  J  H  o'  ij 


Here,  the  induced  dipole  at  force  center  i  is  *j  =  •  iC]0*,  • ,  is 


the  isotropic  atomic  polarizability,  C}01  is  the  total  electro¬ 
static  field  at  the  atomic  site  i  due  to  permanent  charges  qj 
and  induced  dipoles  •j,  *0  is  the  dielectric  permittivity  of 
vacuum,  is  the  electric  field  due  to  fixed  charges  only,  A. . 
and  B..  are  the  repulsion  parameters,  and  C..  is  the  disper¬ 
sion  parameter  for  interaction  between  atoms  i  and  j  that 
have  atom  types  •  and  • .  The  term  D*12/B..r|j*1z,  with  D 
=  5*  10*5  kcal/mol  for  all  pair  interactions,  is  essentially 
zero  at  typical  atomic  separations  but  becomes  the  dominant 
term  at  r^*  1  A  ensuring  that  UR0*r*  is  repulsive  at  dis¬ 
tances  much  smaller  than  the  size  of  an  atom.  The  Thole 
screening33,34  *aT=0.2*  that  smears  induced  dipoles  has  been 
used  in  order  to  prevent  a  polarization  catastrophe  from  oc¬ 
curring.  Finally,  for  cross-term  •heteroatom*  interactions,  the 
modified  Waldman-Hagler  combining  rules34,35  were  used 


A, ,  -  *  A  .A, 


-1/6 


■'ll  -  1  ,  JU  R 3  D 3  ’ 


B.J- 


b;,6 + b;,6 


Cij=  CijCjj. 


Two  FFs  for  TATB  were  developed  in  this  study:  a  polariz¬ 
able  FF  with  no  intramolecular  nonbonded  interactions  and  a 
nonpolarizable  model  that  similarly  does  not  include  any  in¬ 
tramolecular  nonbonded  interactions.  The  necessity  of  ex¬ 
cluding  intramolecular  nonbonded  interactions  is  discussed 
below. 


B.  Force  field  parameterization 

In  this  work  we  followed  the  same  protocol  that  has  been 
applied  previously  to  FF  development  for  a  number  of  mo¬ 
lecular  systems31 34  including  our  recent  work  on  pentaeryth- 
ritol  tetranitrate  *PETN*.32  Here  we  only  briefly  outline  the 
procedure  and  discuss  details  specific  to  TATB.  First,  atomic 
polarizabilities  are  determined  by  fitting  to  the  molecular  po¬ 
larizability  of  gas-phase  molecules  as  determined  from  QC 
calculations.  Second,  partial  atomic  charges  are  fit  to  repro¬ 
duce  QC  results  for  the  electrostatic  potential  on  a  grid  of 
points  around  a  gas-phase  molecule,  while  simultaneously 
providing  a  good  description  of  molecular  dipole  and  quad¬ 
ruple  moments  *see  Eq.  *5*  below*.  Third,  bond  lengths  and 
equilibrium  bending  angles  are  adjusted  to  reproduce  geom¬ 
etries  of  gas-phase  molecules  from  QC  while  bending  force 
constants  are  taken  either  from  previously  developed  FFs  or 
from  fits  to  QC  energies  for  bending-angle  distortions  spe¬ 
cific  to  TATB.  Finally,  parameters  for  torsional  and  out-of¬ 
plane  deformations  are  determined  by  fitting  the  gas-phase 
conformational  energies  obtained  from  QC  calculations  for 
selected  distortions  in  TATB 

Atomic  polarizabilities.  Calculations  of  gas-phase  atomic 
polarizabilities  for  TATB,  performed  at  the  M052X/aug-cc- 
pvDz  level  of  theory  and  basis  set,36  indicate  very  large  an¬ 
isotropy;  specifically,  the  molecular  polarizability  in  the  di¬ 
rection  perpendicular  to  the  molecujar  plane  is  about  a  factor 
of  three  smaller  than  the  in-plane  polarizability  •*. 
=  10.9  A3,  *.=30.6  A3*.  However,  in  the  FF  form  used 
here,  molecular  polarizability  is  determined  as  a  sum  of  iso- 
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TABLE  I.  Force  field  parameters. 


Atom  type 

Atomic 

polarizability 

•  -A3- 

Atomic  charge 
q  *e* 

A  B 

•kcal  mol’1*  -A'1* 

C 

•kcal  A6  morT* 

C""  no2- 

0.70 

0.2839 

107  023.9  3.640S 

$$4.01 

o  nh2- 

0.70 

•  0.202$ 

107  023.9  3  640$ 

$$4.01 

N**  02* 

1.0 

0.S66S 

13  783.9  3.3333 

423.36 

N"  H2* 

1.0 

•  0.  $77$ 

39  091.8  3.31S8 

833.48 

0 

1.0 

•  0.42S2 

15  923.10  3  6446 

239.07 

H 

0.21 

0.3900 

7584.2  5.2846 

8.23 

Bond  length 

Bend  type 

Kb 

•o 

Bond  type 

•constrained,  A* 

Ub^O.SKb-**  V2 

•kcal  mol’ 1  rad2* 

•deg* 

C-C 

1.440 

C-ONIV-C 

144 

118.4 

c-n-o2- 

1.427 

C-ONOp— C 

144 

121.6 

c-n*h2* 

1.322 

C-C-N«H2- 

140 

121.7 

N-0 

1.230 

c-c-no2« 

140 

119  2 

N-H 

1.010 

C-N-0 

140 

116.8 

ON-H 

no 

112.0 

O-N-O 

9$ 

117.3 

H-N-H 

12S 

119.0 

Torsion  type 

Kt,i 

*1.2 

Improper  dihedral  type 

*ID 

Uj  =  •  O.SKT  ,*1  •  cos*i*  •• 

•kcal  mol’1* 

•kcal  mol’ 1# 

Uio=0.SK|q-  2 

•kcal  mol’ 1  rad2* 

c-c-c-c 

0.0 

$.3$ 

C-C-C-N* 

36.  S 

C-C-C-N 

0.0 

20.0 

O-N-0 -C 

89.3 

N-C-C-N 

0.0 

20.0 

h-n-h-c 

2.1 

H-N-C-C 

0.0 

8.88 

O-N-C-C 

0.0 

1.60 

tropic  atomic  polarizabilities,  which  leads  to  isotropic  over- 
allmolecular  polarizability  because  intramolecular  induced- 
dipole/induced-dipole  interactions  are  neglected.  The  latter 
choice  was  made  in  order  to  improve  the  description  of  the 
electrostatic  potential  around  a  molecule  -as  described  be¬ 
low*.  Thus,  while  we  can  match  the  scalar  value,  averaged 
over  all  directions,  of  molecular  polarizability  obtained  from 
QC  it  is  impossible  within  our  approach  to  capture  the  an¬ 
isotropy  of  the  molecular  polarizability.  In  previous  works 
we  assigned  atomic  polarizabilities  such  that  the  average  mo¬ 
lecular  polarizability  is  equal  to  that  obtained  from  QC  cal¬ 
culations.  Initially,  we  followed  the  same  procedure  for 
TATB;  some  atomic  polarizabilities  were  adopted  from  pre¬ 
viously  developed  FF  for  alkanes  and  ethers  while  others 
were  adjusted  to  match  the  average  molecular  polarizability 
•  =23.6  A3.  We  subsequently  determined  that  at  high  pres¬ 
sures  simulations  using  these  polarizabilities  resulted  in  poor 
convergence  of  the  iterative  procedure  used  for  determina¬ 
tion  of  induced  dipole  moments.  We  think  this  was  due  to 
significant  overestimation  of  the  molecular  polarizability  in 
the  direction  perpendicular  to  the  molecular  plane  that  re¬ 
sulted  from  the  use  of  the  average  molecular  polarizability  in 
those  simulations.  This  perpendicular  direction  corresponds 
to  the  stacking  motif  in  TATB  crystal  *Fig.  1*  and  is  the  one 
that  undergoes  the  largest  changes  under  compression. 
Therefore,  the  atomic  polarizabilities  for  carbon  and  nitrogen 
were  reduced  empirically  to  facilitate  reliable  convergence  of 
induced  dipoles  at  higher  pressure  as  well  as  provide  a  better 


agreement  of  molecular  polarizability  with  that  predicted 
from  QC  calculations  for  the  perpendicular  direction,  *Even 
so,  the  results  discussed  in  Sec.  Ill  below  suggest  that  treat¬ 
ment  of  induced  polarization  requires  additional  adjustment 
for  simulations  at  pressures  in  excess  of  6  GPa.*  The  final 
atomic  polarizabilities  are  given  in  Table  I;  they  result  in  an 
isotropic  molecular  polarizability  of  •  =17.46  A3. 

Partial  atomic  charges.  In  order  to  establish  partial 
atomic  charges,  the  electrostatic  potential  was  calculated  at 
the  MP2/aug-cc-pvDz//M052X/aug-cc-pvDz  level  for  a  Car¬ 
tesian  grid  of  •  105  evenly  spaced  points  surrounding  the 
low-energy  gas-phase  conformer  of  TATB.  The  molecular 
dipole  and  quadrupole  *  j  moments  were  calculated  using 
the  same  combinations  of  theory  and  basis  set.  Charge  bond 
increments  were  used  to  calculate  partial  atomic  charges.  The 
value  of  the  partial  charge  q;  positioned  on  atom  i  is  calcu¬ 
lated  as  a  sum  of  all  charge  bond  increments  that  involve 
atom  i  as  shown  in  Eq,  *4*. 


n  Bonds 

q  =  •  *,j.  *4* 

J 


A  set  of  charge  bond  increments  •  *•  was  determined  by  mini¬ 
mizing  the  objective  function 
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>1  J  J 

+  .  *  oc  .  ,  FF.  I. .2®  .5. 


where  ‘j*  and  •  FF  are  the  electrostatic  potential  at  thejth 
grid  point  from  QC  and  the  trial  FF,  ^  and  *'r  are  mo¬ 
lecular  dipole  moments,  and  •  oc  and  •  FF  are  quadrupole 
moments  from  QC  and  the  FF,  respectively.  The  relative 
weights  •',**,•'  for  fitting  the  electrostatic  potential,  di¬ 
pole  moments,  and  quadrupole  moments,  respectively,  were 
set  to  1.0,  0.1,  and  0.05.  Electrostatic  potential  grid  points 
closer  than  1.8  A  to  oxygen,  1.5  A  to  hydrogen,  2.5  A  to 
carbon,  and  2.2  A  to  nitrogen  atoms  were  excluded  from  the 
fit,  as  were  points  further  than  4.0  A  from  the  nearest  atom. 

In  FFs  for  molecular  systems  the  van  der  Waals  and 
electrostatic  interactions  are  typically  excluded  for  closely 
bonded  atoms  -separated  by  one  or  two  bonds,  the  so-called 
1-2  and  1-3  interactions*.  For  atoms  separated  by  three 
bonds  *1-4  interactions*  these  interactions  are  either  in¬ 
cluded  fully  or  scaled  by  some  factor  depending  on  the  FF  or 
•in  some  cases*  chemical  details  of  the  system.  For  nonpo- 
larizable  FFs  the  fitting  of  partial  atomic  charges  using  the 
procedure  described  above  is  independent  of  the  details  of 
exclusion  of  intramolecular  interactions  since  those  interac¬ 
tions  do  not  influence  the  electrostatic  field  around  the  mol¬ 
ecule.  However,  for  polarizable  FFs  the  details  of  intramo¬ 
lecular  interactions  become  important  since  dipole  moments 
induced  on  each  atom,  which  are  determined  by  both  inter- 
and  intramolecular  interactions,  do  influence  the  electrostatic 
field  at  a  given  grid  point.  Initially,  we  attempted  to  fit  the 
partial  atomic  charges  while  allowing  for  intramolecular  in¬ 
teractions  between  atoms  separated  by  three  or  more  bonds. 
However,  we  found  that  the  objective  function  in  Eq.  *5*  was 
uncharacteristically  large  **‘*2*  =  1.8  kcal/mol  per  grid 
point*,  indicating  a  relatively  poor  description  of  the  electro¬ 
static  potential  around  the  molecule.  As  can  be  seen  from  the 
molecular  structure  *Fig.  1*  TATB  has  several  atomic  pairs 
that  are  separated  by  more  than  three  bonds  yet  have  quite 
small  interatomic  separations.  In  this  case  it  is  difficult  to 
define  a  simple  bond-based  exclusion  list  for  intramolecular 
interactions  that  would  provide  a  reasonably  charge-neutral 
interaction  between  individual  atoms  and  the  remaining  parts 
of  the  molecule.  For  example,  if  we  select  any  nitrogen  atom 
and  apply  standard  exclusion  of  1-2  and  1-3  terms  then  its 
interaction  with  carbon  atoms  of  neighboring  amine  *or  ni- 
tro*  groups  will  be  excluded  while  interactions  with  nitrogen, 
hydrogen  *or  oxygen*  atoms  that  are  further  removed  in 
terms  of  connectivity  yet  are  in  very  close  proximity  spa¬ 
tially  will  be  included.  We  believe  that  this  leads  to  an  incor¬ 
rect  distribution  of  induced  atomic  polarization  that  in  turn 
results  in  a  poor  description  of  the  electrostatic  potential.  The 
more  appropriate  way  to  exclude/include  intramolecular  in¬ 
teractions  would  be  to  use  a  group-based  approach  where 
each  group  is  defined  such  that  its  total  charge  is  approxi¬ 
mately  zero.  However,  while  doable,  implementation  of  such 
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an  approach  is  rather  complicated  taking  into  account  the 
molecular  structure  of  TATB.  Therefore,  we  have  decided  to 
exclude  all  intramolecular  interactions  *van  der  Waals  and 
electrostatic*  in  our  FF  model.  In  this  case,  no  dipole  mo¬ 
ments  are  induced  due  to  intramolecular  interactions  and 
hence  description  of  the  electrostatic  potential  around  a  mol¬ 
ecule  depends  only  on  the  arrangement  of  fixed  partial 
atomic  charges.  Charge  fitting  with  this  approach  yielded  a 
much-improved  description  of  the  electrostatic  potential 
.♦••2.=0.88  kcal/mol  per  grid  point*.  The  resulting  partial 
atomic  charges  obtained  from  the  fit  of  bond  increments 
based  on  this  procedure  are  given  in  Table  I, 

Repulsion-dispersion  parameters.  The  C-C  repulsion 
and  dispersion  parameters  were  fit  to  the  experimental  den¬ 
sity  and  •  of  benzene,  those  for  N*H2*-N*H2*  and  H-H 
were  fit  to  to  the  experimental  densities  and  heats  of  vapor¬ 
ization  for  hydrazine  and  methyl  hydrazine,  those  for  0-0 
were  taken  from  the  previous  work  on  PETN,32  and  those  for 
N*02*-N*02*  were  adjusted  to  optimize  the  description  of 
TATB  unit  cell  parameters  and  •  H^.  In  Ref.  37  we  provide 
a  table  containing  density  and  •  Hvap  obtained  using  these 
repulsion-dispersion  parameters  for  benzene,  hydrazine,  and 
methyl  hydrazine. 

The  nitrogen  atoms  in  amine  groups  and  nitro  groups 
were  treated  separately.  Due  to  the  relatively  large  with¬ 
drawal  of  electron  density  from  the  nitro  group  nitrogen 
atom  by  the  bonded  oxygen  atoms,  that  nitrogen  has  a  rela¬ 
tively  large  positive  charge;  thus,  it  is  reasonable  to  expect 
that  it  will  have  effectively  smaller  van  der  Waals  interac¬ 
tions  compared  to  the  amine  nitrogen.  All  repulsion  and  dis¬ 
persion  parameters  are  given  in  Table  I. 

Valence  bonds  and  bends.  Bond  lengths  were  set  to  the 
average  values  of  those  in  an  isolated  TATB  molecule  as 
predicted  at  the  M052X/aug-cc-pvDz  level.  In  our  previous 
simulations32  of  PETN  crystals  we  showed  that  for  hydro¬ 
static  compressions  of  up  to  9  GPa  *the  highest  pressure 
investigated  in  that  study*  constraining  the  covalent  bond 
lengths  had  little  effect  on  the  predicted  properties  of  the 
crystal.  Since  using  fixed  covalent  bond  lengths  allows  use 
of  a  considerably  larger  time  step  compared  to  simulations 
with  flexible  bonds,  in  the  present  work  we  chose  to  use 
fixed  covalent  bond  lengths;  consequently,  bond  stretching 
force  constants  were  not  parameterized.  Equilibrium  bending 
angles  **0*  were  fit  to  reproduce  the  equilibrium  geometry  of 
TATB  predicted  at  M052X/aug-cc-pvDz  level.  The  mean- 
square  deviation  of  TATB  bending  angles  obtained  from  mo¬ 
lecular  mechanics  calculations  with  the  FF  from  the  corre¬ 
sponding  QC  values  was  0.30°.  Most  of  the  bending  force 
constants  were  transferred  from  previous  FFs11-32  without 
change;  exceptions  are  parameters  for  the  •  NH2  group, 
which  were  fit  to  reproduce  MP2/aug-cc-pvDz//M052X/aug- 
cc-pvDz  energies  for  amine  group  bending  and  out-of-plane 
deformation.  Equilibrium  bond  lengths,  bending  angles  and 
three-center  angle-bending  force  constants  are  given  in 
Table  I. 

Out-of-plane  deformation  potentials.  The  H-N-H-C" 
out-of-plane  deformation  was  adjusted  to  reproduce  QC  cal¬ 
culations  at  the  MP2/aug-cc-pvDz//M05-2X/aug-cc-pvDz 
level  for  distortion  of  the  out-of-plane  angle  given  here  as 
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deformation  angle 


FIG.  2.  Relative  energies  for  -a-  the  deformation  energy  for  H-N-H-C' 
out-of-plane  bending,  *b*  the  •  NH?  group  rotation,  and  *c*  the  •  NO*  group 
rotation  as  obtained  from  molecular  mechanics  -MM*  calculations  using  the 
FF  developed  in  the  present  study  and  QC  calculations  at  the  MP2/aug  cc- 
pvDz//M05-2X/aug-cc-pvDz  level.  In  both  QC  and  MM  calculations  the 
♦  NH;  and  ♦  NO;  groups  have  been  constrained  to  planar  geometry. 


the  angle  between  the  C-N  bond  and  the  NH2  plane.  The 
comparison  between  QC  data  and  FF  prediction  for  potential 
energy  as  a  function  of  the  out-of-plane  angle  is  shown  in 
Fig.  2*a\  The  out-of-plane  bending  potential  for  the 
O-N-O-C’  was  taken  to  be  the  same  as  for  the 
O-N-O-N'  from  our  previous  work  on  DMNA.38  Force 
constants  for  all  out-of-plane  deformations  used  in  this  work 
are  given  in  Table  I. 

Dihedral  potentials.  Dihedral  potentials  were  determined 
by  fitting  to  the  relative  energies  for  independent  rotations 
about  the  C-C-N-0  and  C-C-N-H  dihedrals  in  TATB  ob¬ 
tained  from  QC  calculations  at  the  MP2/aug-cc-pvDz//M05- 
2X/aug-cc-pvDz  level.  The  resulting  dihedral  parameters  are 
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given  in  Table  I.  A  comparison  between  QC  and  FF  predic¬ 
tions  for  the  conformational  energies  and  geometries  for 
these  moieties  is  provided  in  Figs.  2-b-  and  2*c*. 

C.  Simulation  protocol 

Our  simulations  of  TATB  crystals  were  performed  using 
the  molecular  simulation  package  Lucretius,39  which  has 
the  capability  to  handle  induced  polarization  effects.  The  sys¬ 
tems  studied  contained  192  TATB  molecules  corresponding 
to  96  unit  cells  *a  4*  4-  6  supercell*.  Periodic  boundary 
conditions  were  used  in  all  simulations.  Covalent  bond 
lengths  were  constrained  using  the  velocity-Verlet  form  of 
the  SHAKE  algorithm.40  The  Ewald  summation  method  was 
used  for  treatment  of  long-range  electrostatic  forces  between 
partial  charges  and  between  partial  charges  and  induced  di¬ 
poles  *k=73  reciprocal  vectors  and  ♦  =0.23  were  used*.  A 
tapering  function 


was  used  to  drive  the  induced-dipole/induced-dipole  interac¬ 
tions  to  2ero  over  a  separation  interval  •  R=RC*  R^*,  where 
Rc  is  a  cutoff  set  to  11  A;  R,apef=lO  A  is  the  distance  at 
which  the  scaling  begins;  and  x  =  R*  Rtaper.  where  R  is  the 
dipole-dipole  separation  for  Raper*  R*  Rc-  This  function  is 
equal  to  unity  *i.e.,  no  scaling*  at  Rtaper  and  smoothly  ap¬ 
proaches  2ero  *i.e„  complete  screening*  at  Rc.  The  derivative 
of  this  function  also  approaches  zero  at  Rc,  therefore  ensur¬ 
ing  that  there  are  no  force  discontinuities  at  the  cutoff  radius. 
Induced  dipoles  were  calculated  via  a  direct  iteration  with  a 
predictor-corrector  method.  A  cutoff  of  11  A  was  used  for  all 
van  der  Waals  interactions  and  the  real  part  of  electrostatic 
interactions  in  the  Ewald  summation. 

The  simulations  were  performed  using  the  combined 
MD-MC  approach  used  in  our  previous  simulations  of  crys¬ 
tals  of  energetic  materials 374  Specifically,  a  500  fs  trajec¬ 
tory  segment  of  isochoric-isothermal  MD  *NVT-MD*  was 
followed  by  20  attempted  changes  in  the  volume  and  shape 
of  the  simulation  cell  based  on  NPT  rigid  molecule  MC 
•NPT-MC*.  The  final  configuration  from  the  kth  NPT-MC 
sequence  and  final  velocities  from  the  kth  NVT-MD  segment 
defined  the  initial  point  in  phase  space  for  k+lst  NVT-MD/ 
NPT-MC  segment.  For  NVT-MD  segments  a  multi -timestep 
reversible  reference  system  propagator  algorithm42  was  em¬ 
ployed;  a  time  step  of  0.5  fs  was  used  for  bonding,  bending, 
dihedral,  and  out-of-plane  deformation  motions,  while  a 
1.0  fs  timestep  was  used  for  all  nonbonded  interactions. 
Starting  with  a  configuration  corresponding  to  the  experi¬ 
mental  structure  at  room  temperature  and  P=0  GPa,  equili¬ 
bration  runs  of  length  200  ps  *8000  NPT-MC  samples*  were 
performed  for  each  thermodynamic  state  considered.  The  ini¬ 
tial  phase  space  point  for  equilibration  at  successively  higher 
pressures  was  taken  from  the  equilibration  run  at  the  next 
lower  pressure  once  it  had  reached  steady-fluctuating  values 
of  the  lattice  parameters;  a  similar  approach  was  used  for 
successive  changes  in  temperature  at  atmospheric  pressure. 
The  MC  step  size  for  a  given  thermodynamic  state  was  ad¬ 
justed  during  the  first  100  ps  of  equilibration  to  yield 
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TABLE  II.  Unit  cell  parameters  and  volume  and  heat  of  sublimation  obtained  from  MD  MC  simulations  and  experiments  at  300  K  and  atmospheric  pressure. 


Experiment 

Simulation, 
polarizable  FF 

Deviation 

Simulation, 
nonpolarizable  FF. 

•polarizable  FF  charges  scaled  by  1.20* 

Deviation 

•%• 

a  *A* 

9.010* 

8946 

•  0.71 

8.984 

*  0.29 

b  *A* 

9.028' 

8.974 

•  0.60 

9.022 

•  0  07 

c  *A* 

6.812* 

6.8S2 

0.S9 

6.86S 

0.78 

•  ♦deg* 

108.S8* 

107.48 

•  1.01 

107.67 

*  0.84 

•  *deg* 

91.82* 

91. 6S 

•  0.18 

90.93 

•  0  97 

*  *deg* 

119.97* 

120.03 

00$ 

120.01 

0.03 

V  -A3* 

442.S2* 

443.4S 

0.21 

449.46 

1.S6 

•  Hsub  ‘kcal/mol* 

43.16 

40.4 

•  6.S 

41.6 

•  3.7 

40.2C 

O.S 

3.S 

'From  Ref.  10. 
“’From  Ref.  23. 
'From  Ref.  24. 


♦  40%-50%  acceptance  probability,  after  which  it  was  held 
constant  for  the  duration  of  equilibration  and  production 
simulation  for  that  state.  Production  runs  consisted  of  1.0  ns 
NVT-MD  *40  000  NPT-MC  observations*  for  each  thermo¬ 
dynamic  state  considered.43 

To  verify  that  the  results  of  our  NVT-MD/NPT-MC 
simulations  are  insensitive  to  particular  choices  of  simulation 
protocol  parameters  we  conducted  four  case  studies  for  one 
of  the  systems.  Specifically,  for  a  nonpolarizable  model  with 
charges  increased  by  20%  relative  to  values  used  in  the  po¬ 
larizable  FF  *see  Sec.  Ill  A*,  simulations  at  300  K  and  atmo¬ 
spheric  pressure  were  performed  using  the  following  combi¬ 
nations  of  simulation  parameters:  *a*  500  fs  NVT-MD,  20 
NPT-MC  attempts  with  *h=0.05  A  *i.e.,  the  original 
scheme  described  above;  here  and  below,  •  h  is  the  maxi¬ 
mum  allowed  change  in  any  of  the  six  nonzero  elements  of 
the  3*  3  upper  triangular  matrix  that  specifies  the  dimen¬ 
sions  of  the  simulation  cell*;  *b*  500  fs  NVT-MD,  20 
NPT-MC  attempts,  *h  =  0.1  A;  *c*  2000  fs  NVT-MD,  80 
NPT-MC  attempts,  •  h=0.05  A;  and  *d*  2000  fs  NVT-MD, 
20  NPT-MC  attempts,  •  h=0.05  A.  By  increasing  •  h  from 
0.05  A  *case  a*  to  0.1  A  *case  b*  the  probability  of  accepting 
a  MC  shape/volume  change  decreased  from  •  48  to  »  23%. 
In  cases  c  and  d,  we  increased  the  length  of  MD  segments 
from  500  fs  to  2  ps  and  either  increased  the  number  of  MC 
attempts  to  80  -thereby  maintaining  the  original  relative  ratio 
of  MD  time  steps  to  MC  trial  moves,  case  c*,  or  kept  the 
number  of  attempted  MC  volume  changes  fixed  at  20  but 
increased  the  total  length  of  NVT-MD  to  4  ns  -case  d*.  In 
each  case  considered,  40  000  MC  volume  changes  were  at¬ 
tempted  and  used  to  form  the  elastic  tensor.  In  Ref.  37  we 
report  the  resulting  elastic  stiffness  coefficients,  as  well  as 
bulk  and  shear  moduli,  obtained  for  each  of  these  four  case 
studies.  Since  the  mechanical  properties  are  defined  in  terms 
of  fluctuations  in  volume  and  shape  of  the  simulation  cell 
they  are  quite  sensitive  to  the  distributions  of  microscopic 
strain  -compared  to  simple  first  moments  used  to  obtain,  for 
example,  the  lattice  parameters  or  simulation  cell  volume*. 
Table  S2  shows  that  for  the  four  cases  studied  19  of  the  21 
elements  of  the  stiffness  tensor  overlap  to  within  two  stan¬ 


dard  deviations  of  the  mean,  while  for  the  isotropic  moduli 
three  out  of  four  properties  yield  agreement  to  within  two 
standard  deviations  of  the  mean. 

Gas-phase  simulations  were  conducted  using  a  Brownian 
dynamics  integrator  for  an  ensemble  of  TATB  molecules 
with  intermolecular  interactions  turned  off  *i.e.,  an  ideal  gas*. 
Results  of  these  simulations  were  used  to  calculate  the  heat 
of  sublimation. 

III.  RESULTS  AND  DISCUSSION 

A.  TATB  unit  cell  parameters  at  ambient  conditions 

Average  unit  cell  parameters  obtained  for  TATB  at  atmo¬ 
spheric  pressure  and  300  K  using  the  polarizable  FF  are  re¬ 
ported  in  Table  II.  where  we  also  provide  experimental  val¬ 
ues  obtained  by  Cady  and  Larson  :  Deviations  for  the  a,  b, 
and  c  unit  cell  parameters  are  less  than  1%,  which  is  com¬ 
parable  to  the  accuracy  of  simulations  using  Gee  et  al.19  and 
TraPPE  *Ref.  22*  nonpolarizable  FFs.  Predicted  unit  cell 
angles,  •  in  particular,  are  in  good  agreement  with  experi¬ 
ment  and  are  somewhat  better  than  those  predicted  from  pre¬ 
vious  simulations.  The  resulting  average  volume  of  the  unit 
cell  is  about  0.2%  larger  than  experimental  value.  Table  II 
also  contains  a  comparison  of  the  heat  of  sublimation  •H* 
obtained  from  our  simulations  to  the  two  experimental  val¬ 
ues.  Unlike  predictions  from  simulations  using  the  FFs  of 
Gee  et  al.19  and  the  TraPPE,  both  of  which  underestimated 
♦  Hsub  by  30%-40%,  our  value  agrees  with  experiment  to 
within  a  few  percent  suggesting  that  our  FF  provides  a  more 
accurate  description  of  the  strength  of  intermolecular  inter¬ 
actions  of  TATB  molecules  in  the  crystal. 

To  investigate  the  importance  of  inclusion  of  induced 
polarization  in  our  simulations  we  also  conducted  simula¬ 
tions  in  which  atomic  polarizabilities  of  all  atoms  were  set  to 
zero.  Since  all  intramolecular  interactions  were  excluded  in 
our  parameterization  of  partial  atomic  charges,  setting  atomic 
polarizabilities  to  zero  does  not  influence  the  description  of 
the  electrostatic  potential  around  a  single  TATB  molecule  in 
the  gas  phase.  All  other  FF  parameters  were  kept  exactly  the 
same  as  in  polarizable  version  of  FF.  Our  simulations  using 
this  nonpolarizable  FF  predicted  cell  dimensions  that  are  off 
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by  +0.54,  +0.73,  and  +0.57%  from  experiments  for  a,  b,  and 
c,  respectively.  Unit  cell  angles  were  106.7°,  92.1°,  and 
120.0°,  and  the  unit  cell  volume  was  456.82,  which  is  3.2% 
greater  than  the  experimental  value  of  Cady  and  Larson.  The 
heat  of  sublimation  obtained  from  simulations  using  this 
nonpolarizable  FF  was  35.9  kcal/mol,  which  is  significantly 
lower  than  the  experimental  value. 

In  principle,  one  can  make  empirical  adjustments  to  the 
FF  to  improve  prediction  from  simulations  using  the  nonpo¬ 
larizable  FF.  In  previous  work  such  empirical  adjustment 
was  implemented  by  uniform  rescaling  of  partial  atomic 
charges,  resulting  in  an  increase  in  the  molecular  dipole  mo¬ 
ment  and  therefore  effectively  compensating  for  the  missing 
♦in  nonpolarizable  FF*  induced  polarization  effects  in  the 
condensed  phase.  Note  that  this  scaling  comes  at  the  expense 
of  a  serious  degradation  in  the  description  of  the  electrostatic 
potential  around  a  single  gas-phase  TATB  molecule  as  well 
as  the  gas-phase  dipole  and  quadrupole  moments.  We  have 
tried  a  similar  approach  for  TATB  crystal,  increasing  the  par¬ 
tial  atomic  charges  by  as  much  as  30%.  We  found  that  simu¬ 
lations  using  the  nonpolarizable  FF  with  partial  charges  in¬ 
creased  by  20%  showed  the  best  ‘among  nonpolarizable 
simulations  with  various  uniform  charge  scalings*  agreement 
with  experiment  as  illustrated  in  Table  II.  As  one  can  see 
from  the  table,  simulations  with  charges  increased  by  20% 
yield  good  reproduction  of  •  but  still  overestimate  the 
unit  cell  volume  by  1.6%.  Further  increase  in  partial  atomic 
charges  did  not  improve  the  description  of  the  unit  cell  vol¬ 
ume  but  continued  to  increase  the  calculated  ♦  H**.  There¬ 
fore,  it  appears  that  improving  prediction  of  both  density  and 
heat  of  sublimation  for  the  nonpolarizable  FF  is  challenging. 
This  further  illustrates  that  inclusion  of  induced  polarization 
is  important  for  capturing  all  important  interactions  in  TATB 
crystal,  which  should  lead  in  turn  to  more  accurate  descrip¬ 
tion  of  structural  and  thermodynamic  properties  of  this  ma¬ 
terial. 


B.  Pressure  dependence  of  TATB  unit  cell  parameters 

MD-MC  simulations  of  the  TATB  crystal  were  con¬ 
ducted  using  polarizable  and  nonpolarizable  *with  20% 
scaled  charges*  FFs  for  pressures  up  to  10  GPa.  In  Fig.  3  we 
show  the  pressure  dependence  of  unit  cell  dimensions.  Re¬ 
sults  obtained  from  simulations  using  both  polarizable  and 
nonpolarizable  FFs  are  compared  to  experimental  data  of  01- 
inger  and  Cady  2  and  Stevens  et  al.  The  latter  are  based 
upon  additional  analysis  by  the  authors  of  Ref.  13  that  did 
not  appear  in  their  original  publication.  For  convenience, 
these  data  are  also  shown  in  tabular  form  in  Table  III.  Figure 
3  reveals  that  despite  scatter  in  experimental  data  there  is 
some  systematic  difference  in  pressure  dependence  of  unit 
cell  dimensions  from  the  two  experimental  data  sets.  Taking 
into  account  this  discrepancy  between  experimental  data, 
predictions  from  simulations  using  both  FFs  are  consistent 
with  experiment.  We  note  that  the  pressure  dependence  of 
the  a,  b,  and  c  unit  cell  parameters  obtained  from  simulations 
using  the  nonpolarizable  FF  is  basically  the  same  as  obtained 
from  the  polarizable  FF. 

Figure  4  compares  the  dependence  of  unit  cell  angles  on 


FIG.  3.  Unit  cell  dimensions  as  a  function  of  pressure  as  obtained  from 
MO-MC  simulations  and  two  experiments. 

pressure  as  obtained  from  simulation  and  experiments.  Note 
that  the  two  sets  of  experimental  results  exhibit  noticeable 
differences;  in  particular,  the  measurements  of  Stevens  et 
al.13  show  noticeably  larger  pressure  dependence  for  •  and  • 
angles  than  in  the  older  measurements.  Simulation  results  are 
in  good  agreement  with  experiment  for  the  angles  •  and  ♦ 
but  systematically  underestimate  ♦.  Finally,  in  Fig.  5  unit 
cell  volume  is  shown  as  function  of  pressure.  Excellent 
agreement  between  both  simulation  results  and  both  sets  of 
experimental  data  is  obtained  over  entire  interval  of  pres¬ 
sures  investigated. 

C.  Temperature  dependence  of  TATB  unit  cel) 
parameters 

Simulations  in  the  temperature  range  between  200  and 
400  K  at  atmospheric  pressure  showed  almost  linear  depen¬ 
dence  of  unit  cell  dimensions  on  temperature.  Simulations 
with  both  the  polarizable  and  nonpolarizable  *20%  scaled 
charges*  FFs  showed  that  thermal  expansion  in  the  direction 
perpendicular  to  TATB  layers  stacking  *along  c  axis*  is  al¬ 
most  an  order  of  magnitude  larger  than  for  two  other  direc¬ 
tions.  For  both  FFs  simulations  predicted  linear  thermal  ex¬ 
pansion  of  1*  10‘4, 1*  I0‘4and9*  10'4  A/K  fora,  b,  and 
c  unit  cell  dimensions,  respectively.  These  thermal  expansion 
coefficients  are  very  similar  to  those  obtained  from  previous 
simulations  using  Gee  et  al. 3  and  TraPPE  *Ref.  22*  FFs. 
Comparison  with  experimental  values  for  linear  thermal  ex¬ 
pansion  coefficients  —  1-  10*4,  2*  10'4,  and  1.7 
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TABLE  III.  Unit  cell  dimensions  extracted  from  experiments  of  Stevens  et  al ,  Ref.  13. 


Pressure 

a 

b 

C 

. 

. 

. 

•GPa* 

•A- 

•A- 

-A- 

•deg* 

•deg- 

•deg* 

0.5600 

8.9671 

9.0821 

6.6242 

110.6190 

93.0440 

118.9250 

0.8500 

8.9372 

90472 

6  5919 

110  4758 

93.2130 

1189100 

1.2000 

8.9201 

9.0250 

6.5571 

110.4850 

93.2510 

118.9660 

1.6600 

8.8846 

9  0001 

6  5199 

1106200 

93.4120 

118.8670 

1.9500 

8.8749 

8.9609 

6  4976 

110.7200 

93.3560 

119.1630 

2.4200 

8.8525 

89424 

6.4607 

110  8767 

93.5040 

118.9640 

3.2900 

8  8087 

8  8885 

6.3963 

111  0133 

93.6910 

119.0320 

4.9500 

8.7373 

8.8547 

6.2898 

111.3570 

94.2640 

118.4730 

8.1800 

8.6510 

8.8133 

6.1145 

112.5060 

95.5259 

117.5540 

8.5300 

8.6247 

8.7889 

6.0885 

112.4983 

95.5370 

117.6220 

9.5300 

8.6027 

8.7814 

6  0461 

112.8150 

95.8440 

117.3280 

10.1700 

8.5716 

8.7608 

6.0225 

112.9080 

96.2880 

116.9200 

11.6400 

8.5433 

8.7458 

5.9713 

113.4380 

96.5210 

116.8300 

13.2200 

8.5073 

8.7629 

5.9309 

114.2140 

96.5930 

116.6340 

•  10' 3  A/K  for  a,  b,  and  c  parameters,  respectively  shows 
that  while  for  the  a  dimension  simulations  somewhat  over¬ 
estimate  the  thermal  expansion  coefficient,  the  opposite  trend 
is  observed  for  the  b  and  c  directions.  Similarly,  simulations 
•both  FFs*  yield  a  volumetric  thermal  expansion  coefficient 
of  0.07  A3/K  which  is  smaller  than  the  experimental  value 
of  0.14  A3/K  It  is  interesting  to  point  out  that  indepen¬ 
dently  of  which  FF  is  used  *Gee  et  al.;19  TraPPE;  this  work, 
polarizable,  or  nonpolarizable  with  increased  charges*,  the 
predicted  thermal  expansion  coefficients  are  very  similar,  in¬ 
dicating  a  surprising  insensitivity  of  thermal  expansion  to  the 
details  of  FF  parameters. 

D.  TATB  elastic  coefficients 

Complete  isothermal  second-order  elastic  tensors  and 
derived  isotropic  bulk  and  shear  moduli  *K  and  G,  respec¬ 
tively*  were  determined  for  TATB  by  applying  the 
Parrinello-Rahman  strain-strain  fluctuation  formula4  to  our 
simulation  data  Results  for  both  FF  models  'polarizable  and 
nonpolarizable  with  scaled  charges*  were  obtained  for  six 
evenly  spaced  pressures  between  0  and  10  GPa  on  the  300  K 
isotherm  and  for  five  evenly  spaced  temperatures  between 
200  and  400  K  at  atmospheric  pressure.  In  the  Parrinello- 


Rahman  approach  the  elastic  compliance  tensor  Sy  is  con¬ 
structed  from  observations  of  simulation  cell  size  and  shape 
obtained  from  an  NIPT  simulation;  the  elastic  stiffness  tensor 
Cy  is  simply  the  matrix  inverse  of  the  compliance  tensor. 
Details  of  our  implementation  of  the  approach  have  been 
published  previously.41  In  the  triclinic  TATB  crystal  all  21 
unique  elements  of  the  elastic  tensor  are,  in  general,  nonzero, 
and  are  reported  here.  To  the  best  of  our  knowledge  there  are 
no  previous  predictions  or  measurements  available  for  com¬ 
parison,  with  the  exception  of  bulk  modulus  values  based  on 
EOS  fits  to  isothermal  compression  data. 

Results  for  300  K  and  1  atm.  Results  for  the  isothermal 
elastic  stiffness  tensor  obtained  at  300  K  and  atmospheric 
pressure  are  reported  in  Table  IV.  Results  for  the  polarizable 
and  nonpolarizable  FFs  are  comparable,  as  is  shown  graphi¬ 
cally  for  the  stiffness  tensor  in  Fig.  6.  The  most  obvious 
difference  seen  there  is  that  Cn  and  C22  for  the  two  FFs 
differ  by  more  than  one  estimated  standard  deviation  of  the 
mean  *1  •  *,  with  larger  values  obtained  for  the  polarizable 
FF.  *ln  fact,  eight  of  the  nine  "orthotropic"  elements  of  the 
stiffness  tensor  are  larger  for  the  polarizable  FF  than  for  the 
nonpolarizable  one,  but  five  of  the  nine  pairs  overlap  within 
1  • ,  and  all  but  Cn  and  C22  overlap  within  2*  .*  The  hydro¬ 
static  compression  and  thermal  expansion  of  TATB  are 
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FIG.  4.  Unit  cell  angles  as  a  function  of  pressure  as  obtained  from  MO  MC  FIG.  5.  Unit  cell  volume  as  a  function  of  pressure  as  obtained  from  MO  MC 
simulations  and  two  experiments.  simulations  and  two  experiments. 
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TABLE  IV.  Elastic  stiffness  coefficients  and  derived  isotropic  moduli  calculated  at  T- 300  K,  P-0  GPa. 
•Units  are  GPa.- 


Polarizable  FF 

Nonpolarizable  FF 

Polarizable  FF 

Nonpolarizable  FF 

C„ 

6S.7*  O  S* 

S7.7-  OS 

Ci6 

1.0-  1.0 

0.0-  1.0 

C22 

62.0*  1.0 

S8.0-  1.0 

C24 

0.6-  0.2 

0.6-  02 

C33 

18.3-  OS 

17.0-  O.S 

C2i 

-0.5-  0.2 

-O.S-  0.3 

1.4-  0.3 

1.0-  0.3 

C?6 

1.0-  1.0 

2.0-  08 

C$5 

0.68-  0.06 

06-  0.2 

C34 

0.2-  0.3 

-0.1-  0.4 

Cm 

21  6-  07 

20.3-  0.8 

C35 

-0.4-  0  1 

-0.3-  0  2 

Ciz 

18.5-  0.5 

16.2-  0.7 

C36 

-0.4-  0.7 

-1.0-  O.S 

C13 

4.0-  2.0 

3.2-  0.5 

C45 

0.1-  0.2 

0.01-  0.04 

C23 

SO-  1.0 

S.7-  0.6 

C46 

0.3-  0.2 

-O.S-  0.1 

C14 

-0.2-  0.3 

0.1-  0.1 

C56 

0.4-  0.1 

0.1-  01 

C15 

-1.0-  0.1 

-0.9-  0.2 

Kr**6 

14.0-  1.0 

13.2-  O.S 

Gfieuss1’ 

2.0-  0  2 

1.8-  0.1 

Kvo*gtC 

22.3-  0.3 

20.3-  0.2 

Cyolgt 

12.6-  0.2 

11. S-  0.2 

"Reported  values  were  obtained  from  the  overall  1  ns  simulation  -40  000  MC  observations-.  Uncertainties 
correspond  to  one  standard  deviation  of  the  mean  obtained  from  four  contiguous  2S0  ps  -10  000  MC  observa¬ 
tion-  subsamples. 

^Obtained  from  the  compliance  tensor  -see  Ref.  37-.  Uncertainties  were  obtained  by  application  of  standard 
error  propagation  methods. 

Obtained  from  the  stiffness  tensor.  Uncertainties  were  obtained  by  application  of  standard  error  propagation 
methods 


highly  anisotropic  and  this  is  clearly  reflected  in  the  indi¬ 
vidual  elastic  tensor  elements.  For  example,  for  the  polariz¬ 
able  FF  we  predict  values  Cn  =65.7  GPa,  C22=62  GPa,  and 
C 33=18. 3  GPa,  These  elastic  coefficients  correspond  nomi¬ 
nally  to  deformations  along  the  a,  b,  and  c  crystal  axes, 
respectively,  and  are  consistent  with  the  TATB  crystal  pack¬ 
ing  motif  in  which  layers  of  essentially  planar  molecules, 
held  together  by  hydrogen-bonds  in  the  a-b  plane,  stack  into 
layers  along  the  c  axis  that  are  held  together  by  weak  disper¬ 
sion  interactions.  Similar  remarks  apply  to  the  shear  elastic 
coefficients,  C44  =  1.4  GPa,  C55=0.68  GPa,  and  C66 
=  21.6  GPa  •polarizable  FF*  corresponding  roughly  to  defor¬ 
mations  of  lattice  angles  *,  -,  and  *,  respectively.  The  angle 


FIG.  6.  Comparison  between  elastic  stiffness  tensor  elements  for  TATB 
crystal  computed  at  -T  =  300  K,  P  =  0  GPa-  using  the  polarizable  -ordinate- 
and  nonpolarizable  -abscissa-  FF  models.  Uncertainties  represent  one  esti¬ 
mated  standard  deviation  of  the  mean.  The  diagonal  line  is  included  as  a 
guide  for  the  eyes;  points  falling  on  the  line  correspond  to  equal  values  for 
both  FF  models. 


*  is  defined  by  the  a  and  b  crystal  axes;  the  energy  required 
to  impose  a  shear  strain  in  the  hydrogen-bonded  plane  that 
contains  them  is  much  larger  than  those  required  to  impose 
shear  strains  in  the  planes  whose  definitions  involve  the 
c-axis.  Hence,  one  would  expect  in  agree¬ 

ment  with  the  calculated  result.  Analogous  considerations 
can  be  used  to  explain  the  relative  magnitudes  for  the  dila- 
tional  coefficients,  C12-  *C13,C23*. 

Isotropic  bulk  and  shear  moduli  were  calculated  using 
standard  expressions  for  the  Vfaigt  and  Reuss  bounds,  corre¬ 
sponding  to  conditions  of  uniform  strain  and  uniform  stress, 
respectively.  These  quantities  are  expressed  in  terms  of  the 
stiffness  and  compliance  tensors,  respectively,  but  as  noted 
above  these  are  simply  matrix  inverses  of  one  another.  The 
Reuss  bound  should  always  be  less  than  or  equal  to  the  Vbigt 
bound  and  is  the  one  that  should  be  compared  to  bulk  modu¬ 
lus  values  obtained  from  EOS  fits  to  hydrostatic  compression 
data.  Values  of  the  calculated  isotropic  moduli  are  included 
in  Table  IV.  The  results  for  the  polarizable  and  nonpolariz¬ 
able  FFs  are  quite  similar.  For  the  polarizable  FF  we  ob¬ 
tained  Kg^55=l4  GPa,  Gf£^s  =  2.0  GPa.  and  for  the  nonpo¬ 
larizable  one  K^^'=  13.2  GPa  and  G^^*=1.8  GPa.  These 
numbers  overlap  to  within  the  estimated  uncertainties.  The 
Voigt  average  bulk  and  shear  moduli  are  factors  of  16  and 
6.3  times  larger  than  the  corresponding  Reuss  average  val¬ 
ues.  These  large  differences  for  the  two  bounds  are  indicative 
of  the  large  anisotropy  of  the  crystal. 

The  atmospheric  pressure  isothermal  bulk  modulus  can 
also  be  determined  from  an  EOS  fit  to  unit  cell  volume  ver¬ 
sus  pressure  data  along  the  300  K  isotherm.  For  this  purpose 
we  fit  V»P»  at  300  K  with  a  third-order  Birch-Murnaghan 
EOS45'46  performed  in  the  P-V/V0  plane  with  uniform 
weighting  of  the  simulated  data  points.  In  Table  V  we  pro- 
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TABLE  V.  Isothermal  bulk  modulus  KBM  and  initial  pressure  derivative 
Kp,  M  obtained  from  fits  to  third-order  Birch-Murnaghan  EOS. 


Simulation 

Simulation 

Experiment 

Experiment 

polarizable 

nonpolarizable 

Stevens  et  al. 

Olinger  et  al. 

FF 

FF 

•Ref.  13*1 

•Ref.  12-b 

KB  M  -GPA- 

14.8 

11.7 

13.6*  0.6 

16.7 

Kfi.M 

8.7 

11.9 
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'’Parameters  originally  reported  in  Ref.  13. 


vide  the  corresponding  isothermal  bulk  moduli  and  initial 
pressure  derivatives  K»-dK/dP.  The  predicted  values  are 
K|5j,=14.8  GPa  and  K™jrl  =  11.7  GPa  for  the  initial  bulk 
moduli  for  the  polarizable  and  nonpolarizable  FFs,  respec¬ 
tively,  and  =8.7  and  K*S5pol=11.9  for  the  correspond¬ 
ing  initial  pressure  derivatives.  The  EOS-based  prediction  of 
bulk  modulus  for  the  polarizable  FF  falls  within  the  •  1  • 
standard  error  envelope  for  the  Reuss  average  value  obtained 
directly  from  the  elastic  compliance  tensor,  whereas  for  the 
nonpolarizable  FF  the  EOS-based  value  is  -  3*  below  the 
Reuss  average  one. 

We  also  include  in  Table  V  the  bulk  modulus  values 
extracted  from  the  experimental  room-temperature  isotherms 
of  Stevens  et  al.13  and  Olinger  et  al.12  by  fitting  to  the  third- 
order  Birch-Murnaghan  EOS.  In  their  paper  Stevens  et  al.13 
included  fits  for  three  different  EOSs— Murnaghan,  Vinet, 
and  third-order  Birch-Murnaghan— and  also  results  reported 
by  or  derived  from  earlier  theoretical  studies  due  to  Pastine 
and  Bernecker  6  and  Byrd  et  al.;17  we  have  not  included  the 
two  theoretical  values  in  the  present  study.  Stevens  et  al.13 
obtained  values  Ke_M=13.6  GPa  for  the  fit  to  their  own  data 
•restricted  to  pressures  below  8  GPa*  and  KB.M=16.7  GPa 
for  their  fit  to  the  isotherm  of  Olinger  et  al.12  which  extended 
up  to  P  =  7.47  GPa. 

Temperature  dependence  of  stiffness  coefficients.  In 
Fig.  7  we  show  for  both  FFs  the  predicted  temperature  de¬ 
pendence  of  the  nine  "orthotropic”  stiffness  coefficients  for 
compression  Cn,  C22,  C33  -panel  a*,  shear  C44,  C55, 
•panel  b\  and  dilation  C12,  C13,  and  C23  -panel  c*.  The  re¬ 
sults  for  the  two  FFs  are  similar  and  qualitatively  reasonable. 
It  is  seen  in  panel  a  that  the  compressive  elements  all  un¬ 
dergo  thermal  softening,  with  the  larger  coefficients  Cu  and 
C22  undergoing  larger  net  decreases  than  C33  over  the  200  K 
temperature  interval  studied.  The  shear  elements  C44  and  C55 
are  nearly  independent  of  temperature  for  both  FFs,  whereas 
Cgg  decreases  by  several  GPa;  starting  from  values  that  differ 
by  •  3  GPa  at  200  K,  the  decrease  in  is  greater  for  the 
polarizable  FF  than  for  the  nonpolarizable  one  such  that  the 
two  values  are  comparable  at  400  K.  All  three  dilational 
stiffness  coefficients  C12,  C13,  and  C23  decrease  with  increas¬ 
ing  temperature;  we  find  decreases  of  •  one  GPa  for  C13  and 
C23,  while  C12  decreases  by  5  GPa  over  the  same  tempera¬ 
ture  interval. 

Pressure  dependence  of  the  stiffness  coefficients.  The 
pressure  dependence  of  selected  stiffness  coefficients  at 
300  K  is  shown  for  both  FFs  in  Fig.  8.  The  changes  with 
pressure  between  zero  and  10  GPa  are  significantly  greater 
than  was  the  case  for  temperature  dependence  over  the  200  K 


FIG  7.  Temperature  dependence  of  selected  stiffness  coefficients  for  TATB 
crystal  at  P^0  GPa  for  the  polarizable  and  nonpolarizable  FFs  -dashed  and 
solid  connecting  lines,  respectively*.  Panel  a:  compression  coefficients  C„, 
C22.  and  Cjj.  Panel  b:  shear  coefficients  C<4.  CS5.  and  C«.  Panel  c:  dilational 
coefficients  Cu,  Cij,  and  C>3.  Uncertainties  represent  one  estimated  stan 
dard  deviation  of  the  mean. 

interval  studied.  With  an  important  exception  noted  below, 
the  results  are  in  accord  with  expectations  based  on  the  dis¬ 
cussion  of  the  general  form  of  the  elastic  tensor  for  TATB 
crystal.  For  example,  C13  and  C22  both  increase  from  •  60  to 

•  250  GPa  under  compression,  whereas  C33  increases  from 

•  20  GPa  up  to  only  •  100  GPa;  and  the  increase  in  is 
much  larger  than  the  increase  in  C44  and  C5S  which  remain 
below  5  GPa  for  all  pressures.  A  large  anisotropy  is  also 
noted  for  the  dilational  coefficients,  with  a  much  larger  in¬ 
crease  in  C12  than  for  C,3  or  C23, 

We  note  that  the  polarizable  FF  yields  seemingly  anoma¬ 
lous  results  for  pressures  above  6  GPa;  seven  of  the  nine 
stiffness  coefficients  shown  in  Fig.  8  soften  considerably  for 
P*  8  GPa,  in  some  cases  to  values  less  than  those  obtained 
at  6  GPa.  For  instance,  at  P=6  GPa,  C13-  200  GPa  for 
both  FFs;  at  P  =  10  GPa,  Cti  for  the  nonpolarizable  FF  has 
increased  to  •  250  GPa  -following  a  smooth,  monotonic 
trend  curve-  while  for  the  polarizable  one  Cn  has  remained 
approximately  constant  over  that  4  GPa  pressure  interval,  -In 
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FIG.  8.  As  in  Fig.  7  except  that  the  pressure  dependence  of  selected  stiffness 
coefficients  is  shown  for  300  K. 


some  cases  the  change  in  behavior  occurs  between  6  and 
8  GPa  while  in  others  it  occurs  between  8  and  10  GPa.- 
Since  there  is  no  evidence  for  a  phase  transition  in  the  lattice 
parameters  and  unit  cell  volumes  on  the  300  K  isotherm 
shown  in  Figs.  3-5,  this  behavior  for  the  stiffness  coeffi¬ 
cients  seems  to  be  unphysical.  We  note  that  in  previous  work 
for  PETN  in  which  a  similar  polarizable  FF  development 
protocol  was  applied,  a  structural  phase  transition  was  ob¬ 
served  at  P-  7  GPa  when  using  the  polarizable  FF  devel¬ 
oped  there;  whereas  there  was  no  evidence  for  one  based  on 
the  nonpolarizable  FF  model.  At  such  high  pressures  inter- 
molecular  separations  between  heavy  atoms  -carbons,  nitro¬ 
gens,  and  oxygens-  are  often  less  than  2.7  A.  It  turns  out  that 
for  such  short  separations  the  selected  Thole  screening  func¬ 
tional  with  aT=0.2  is  not  strong  enough  to  sufficiently  damp 
interactions  between  induced  dipoles,  which  leads  to  unreal¬ 
istically  weaker  repulsive  interactions  between  these  atoms 
and  therefore  results  in  softer  stiffness  coefficients  obtained 
from  simulations.  We  found  that  reduction  in  aj  to  0.1,  in 
principle,  eliminates  this  artifact  and  simulations  with  polar¬ 
izable  FF  predict  elastic  constants  that  follow  a  trend  similar 
to  those  obtained  using  nonpolarizable  FF  for  pressures 
above  6  GPa.  However,  while  this  empirical  adjustment  ap- 


FIG.  9.  Temperature  dependence  at  P  =  0  GPa  -panel  a-  and  pressure  de¬ 
pendence  at  T * 300  K  -panel  b-  of  the  Reuss  average  isotropic  moduli  for 
TATB  crystal  for  the  polarizable  and  nonpolarizable  FFs  -dashed  and  solid 
connecting  lines,  respectively-.  Uncertainties  represent  one  estimated  stan¬ 
dard  deviation  of  the  mean. 

pears  to  provide  reasonable  predictions,  aT  =  0.1  in  the  se¬ 
lected  Thole  screening  functional  form  results  in  an  unreal¬ 
istically  large  screening  length  for  the  induced-dipole/ 
induced-dipole  interaction.  We  believe  that  perhaps  use  of 
other  functional  forms  with  more  realistic  values  for  the 
screening  length  '  or  implementation  of  Gaussian  distributed 
♦instead  of  point-  charges  and  dipoles  would  likely  resolve 
this  apparent  roadblock  for  reliable  use  of  polarizable  FF  for 
simulations  at  high  pressures.  We  recommend  that  applica¬ 
tions  of  the  polarizable  FF  with  current  treatment  of  induced 
polarization  be  restricted  to  pressures  of  6  GPa  or  less. 

Temperature  and  pressure  dependence  of  isotropic  bulk 
and  shear  moduli.  Finally,  we  present  in  Fig.  9  the  tempera¬ 
ture  dependence  -panel  a-  and  pressure  dependence  -panel  b- 
of  the  Reuss  average  isotropic  moduli  for  both  FFs.  With  the 
exception  noted  above  for  the  polarizable  FF  at  high  pres¬ 
sure,  the  results  for  the  two  FFs  are  similar.  The  bulk  modu¬ 
lus  for  both  FFs  exhibits  modest  decreases  with  temperature 
•panel  a-.  The  shear  modulus  for  the  nonpolarizable  FF  de¬ 
creases  from  ♦  2.4  to  •  1  GPa  between  200  and  400  K, 
whereas  for  the  polarizable  FF  the  shear  modulus  is  essen¬ 
tially  independent  of  temperature  -panel  a-.  Under  compres¬ 
sion,  the  bulk  modulus  for  the  nonpolarizable  FF  increases 
from  ♦  14  to  •  82  GPa,  whereas  the  shear  modulus  in¬ 
creases  only  from  •  2  to  •  4  GPa  -panel  b-.  The  behavior  of 
the  stiffness  tensor  noted  for  the  polarizable  FF  in  Fig.  8 
results  in  increased  negative  curvature  in  the  bulk  modulus 
for  P-  6  GPa  but  little  if  any  discernable  effect  on  the  shear 
modulus  compared  to  the  result  for  the  nonpolarizable  FF. 

tV.  CONCLUSIONS 

Atomistic  simulations  of  TATB  crystal  have  been  con¬ 
ducted  using  polarizable  and  nonpolarizable  FFs.  We  found 
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that  inclusion  of  induced  polarization  effects  is  essential  in 
providing  simultaneously  accurate  descriptions  of  crystal 
structure  and  heat  of  sublimation  of  TATB  crystal.  However, 
it  appears  that  these  effects  can  be  effectively  taken  into 
account  in  the  nonpolarizable  FF  by  increasing  partial  atomic 
charges  by  20%  relative  to  values  in  the  polarizable  form. 
The  latter  approach  is  about  a  factor  of  two  more  computa¬ 
tionally  efficient  at  normal  pressure  and  a  factor  of  2.6  more 
efficient  at  high  pressure  -10  GPa*  compared  to  the  polariz¬ 
able  model.  Predicted  pressure  dependencies  of  unit  cell  di¬ 
mensions  are  in  very  good  agreement  with  experimental  data 
in  the  entire  pressure  range  up  to  10  GPa.  Heat  of  sublima¬ 
tion  and  temperature  dependence  of  unit  cell  dimensions  ob¬ 
tained  from  simulations  were  also  found  to  be  in  good  agree¬ 
ment  with  experiment.  This  agreement  with  experimental 
data  provides  us  confidence  that  simulation  predicted  elastic 
coefficients  of  TATB  crystal  and  reported  here  for  the  first 
time  are  reliable.  Analysis  of  stiffness  coefficients  showed 
larger  anisotropy  of  elastic,  shear,  and  dilation  components 
consistent  with  the  sheet  stacking  motif  of  TATB  crystal. 
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